/**
 * @file mt2d.cpp
 * @author Yi Zhang (zhangyiss@icloud.com)
 * @brief 2-D Forward modeling of magnetotelluric data at a given frequency
 * @version 1.0
 * @date 2021-12-06
 * 
 * @copyright Copyright (c) 2021
 * 
 * Reference: Key, K., & Weiss, C. (2006). Adaptive finite-element modeling using unstructured grids: 
 * The 2D magnetotelluric example. Geophysics, 71(6), G291–G299. https://doi.org/10.1190/1.2348091
 * 
 * +z downward, z=0 the Air-Earth interface
 * time dependence exp(-iwt)
 * 
 * TE-mode: E = (Ex,0 ,0 ), H = (0, Hy,Hz)
 * TM-mode: E = (0 ,Ey,Ez), H = (Hx,0 ,0 )
 *          x (Ex for TE mode and Hx for TM mode)
 *        /
 *      /
 *    /    air (conductivity = 1e-10, epsilon = 1.0, mu = 1.0, bottom depth = 0)
 *  /
 * O---------------------> y
 * |   ...
 * |---------------------
 * |   layer bottom (conductivity = 0.1, epsilon = 1.0, mu = 1.0, bottom depth = 1e+20)
 * |
 * z
 * 
 */

/*The followings bypass a bug that VsCode's IntelliSense reports incorrect errors when using eigen3 library*/
#if __INTELLISENSE__
#undef __ARM_NEON
#undef __ARM_NEON__
#endif

// Provide FEM shape functions and the Gauss Legendre Numerical Integral (GLNI) functions
#include "../utilities/glni.h"
#include "../utilities/nfem_sf.h"
// Provide primary filed and boundary conditions
#include "../utilities/layerEarth_MT.h"
// Provide array structures and file io
#include "gctl/core.h"
#include "gctl/io.h"
// Provide sparse matrix and solver
#include "Eigen/Dense"
#include "Eigen/Sparse"


#define pi 3.1415926535897932384626433832795
#define _mu0 (4e-7*pi)
#define _eps0 8.8541878176e-12

typedef std::complex<double> complex_d;
typedef Eigen::SparseMatrix<complex_d> spmat_cd; // 注意Eigen默认的稀疏矩阵排序为列优先
typedef Eigen::Triplet<complex_d> triplt_cd;
typedef Eigen::VectorXcd vector_cd;

static const complex_d _one = complex_d(1.0, 0.0);
static const complex_d _zero = complex_d(0.0, 0.0);
static const complex_d _1i  = complex_d(0.0, 1.0);

enum MT_Mode {TE, TM, Modeling};

struct conductive_region
{
    int tag;
    double cndt, mu, epsl; // cndt (conductivity) epsl (epsilon)
    double ref_cndt; // background conductivity
};

struct mt_station
{
    gctl::vertex2dc *vert;
    double app_resist; // apparent resistivity
    double phase; // phase angle
    std::vector<gctl::triangle2d*> host_tris;
};

class mt2d
{
public:
    mt2d(){}
    virtual ~mt2d(){}
    void run(int argc, char *argv[]); // 执行计算流程

protected:
    void read_config(std::string filename); // 读入配置文件
    void read_model(std::string filename); //读入Triangle模型
    void write_model(std::string filename, bool debug = true); // 写入Gmsh模型文件
    void write_observation(std::string filename);

    void get_element_conductivity(int tag_num, double &conduct, bool ref_cndt = false);
    void assemble_linear_system();
    void apply_boundary_condition();
    void solve();
    void post_process();

    // 偏微分方程的左端项
    static double _Aij_part1(double x, double y, double x1, double x2, double x3, 
		double y1, double y2, double y3, void *att, void *instance) // 辅助函数 类的成员函数的外部接口
	{
		return reinterpret_cast<mt2d*>(instance)->Aij_part1(x, y, x1, x2, x3, y1, y2, y3, att);
	}

    static double _Aij_part2(double x, double y, double x1, double x2, double x3, 
		double y1, double y2, double y3, void *att, void *instance)
	{
		return reinterpret_cast<mt2d*>(instance)->Aij_part2(x, y, x1, x2, x3, y1, y2, y3, att);
	}

	double Aij_part1(double x, double y, double x1, double x2, double x3, 
		double y1, double y2, double y3, void *att) // 偏微分方程的左端项1
	{
        unsigned int *uv_id = static_cast<unsigned int*>(att);
        double u = lsf_.triangle(x, y, x1, x2, x3, y1, y2, y3, uv_id[0], nfem_sf::Dx)*lsf_.triangle(x, y, x1, x2, x3, y1, y2, y3, uv_id[1], nfem_sf::Dx);
        double v = lsf_.triangle(x, y, x1, x2, x3, y1, y2, y3, uv_id[0], nfem_sf::Dy)*lsf_.triangle(x, y, x1, x2, x3, y1, y2, y3, uv_id[1], nfem_sf::Dy);
        return (u + v);
	}

    double Aij_part2(double x, double y, double x1, double x2, double x3, 
		double y1, double y2, double y3, void *att) // 偏微分方程的左端项2
	{
		unsigned int *uv_id = static_cast<unsigned int*>(att);
        double u = lsf_.triangle(x, y, x1, x2, x3, y1, y2, y3, uv_id[0], nfem_sf::Value);
        double v = lsf_.triangle(x, y, x1, x2, x3, y1, y2, y3, uv_id[1], nfem_sf::Value);
        return u*v;
	}

    // 偏微分方程的右端项
    static std::complex<double> _Bi(double x, double y, double x1, double x2, double x3, 
		double y1, double y2, double y3, void *att, void *instance)
	{
		return reinterpret_cast<mt2d*>(instance)->Bi(x, y, x1, x2, x3, y1, y2, y3, att);
	}

    std::complex<double> Bi(double x, double y, double x1, double x2, double x3, 
		double y1, double y2, double y3, void *att)
	{
        unsigned int *u_id = static_cast<unsigned int*>(att);
        le_.calculate_EH(-1.0*y, x);
        // 高斯积分是否可以直接应用在复数上还需要验证 但对于变化不大的函数问题不大
        if (work_mode_ == TE) return lsf_.triangle(x, y, x1, x2, x3, y1, y2, y3, *u_id, nfem_sf::Value)*le_.get_Ex();
        else return -1.0*lsf_.triangle(x, y, x1, x2, x3, y1, y2, y3, *u_id, nfem_sf::Dy)*le_.get_Ey(); // 注意纵轴向上与层状模型坐标轴方向相反
	}

    // 以下为类的变量
    unsigned int elem_num_; // 有限元个数（分割的单元数）
    unsigned int node_num_; // 节点个数
    gctl::array<int> node_markers_; // 边界顶点的标签值
    gctl::array<int> elem_markers_; // 元素的标签值
    gctl::array<gctl::vertex2dc> mesh_nodes_; // 网格顶点
    gctl::array<gctl::triangle2d> mesh_elems_; // 网格元素
    gctl::gmshio meshio_; // 网格文件IO

    glni::triangle g2d_tri_; // 二维三角形区域高斯勒让德积分对象
    nfem_sf::linear lsf_; // 线性边矢量形函数对象
    spmat_cd A;
    vector_cd x, B;

    std::string tri_file_; // 网格文件名
    double omega_; // 电磁波角频率
    double theta_in_; // 主场入射角
    int station_marker_; // 观测点标签值
    int bound_marker_; // 边界顶点标签值
    int air_marker_; // 空气层标签值
    std::vector<conductive_layer> primary_model_; // 主场模型 (导电率 相对磁导率 相对介电常数 底深度)
    gctl::array<conductive_region> secondary_model_; // 磁场模型 (标签值 导电率 参考导电率 相对磁导率 相对介电常数)

    int obs_num_;
    std::vector<mt_station> obs_site_; // 观测站

    layerEarthMT le_;

    MT_Mode work_mode_;
};

void mt2d::read_config(std::string filename)
{
    std::clog << "Reading configuration files ... ";
    std::ifstream infile;
    gctl::open_infile(infile, filename, ".config");

    double frequency;
    std::string w_mode;
    infile >> tri_file_ >> w_mode >> frequency >> theta_in_ 
        >> station_marker_ >> bound_marker_ >> air_marker_;

    if (w_mode == "TE") work_mode_ = TE;
    else if (w_mode == "TM") work_mode_ = TM;
    else if (w_mode == "Modeling") work_mode_ = Modeling;
    else throw gctl::runtime_error("Invalid MT surveying mode.");

    omega_ = 2.0*pi*frequency;

    int tmp_size;
    infile >> tmp_size;
    if (tmp_size < 2)
    {
        throw gctl::runtime_error("Invalid primary model sizes: "+std::to_string(tmp_size));
    }

    // Set primary model
    primary_model_.resize(tmp_size);
    for (int i = 0; i < tmp_size; ++i)
    {
        infile >> primary_model_[i].cndt >> primary_model_[i].mu 
            >> primary_model_[i].epsl >> primary_model_[i].btm_depth;
    }

    le_.set_layers(primary_model_);
    if (work_mode_ == TM) le_.set_to_TM_mode();
    le_.set_frequency(frequency);
    le_.initialize();

    infile >> tmp_size;
    if (tmp_size <= 0)
    {
        throw gctl::runtime_error("Invalid secondary model sizes: "+std::to_string(tmp_size));
    }

    // Get secondary model
    secondary_model_.resize(tmp_size);
    for (int i = 0; i < tmp_size; ++i)
    {
        infile >> secondary_model_[i].tag >> secondary_model_[i].cndt >> secondary_model_[i].ref_cndt 
            >> secondary_model_[i].mu >> secondary_model_[i].epsl;
    }

    infile.close();
    std::clog << "done\n";
    return;
}

// Read Triangle's output files
void mt2d::read_model(std::string filename)
{
    std::clog << "Reading Triangle files ... ";
    gctl::matrix<double> tmp_markers;
    gctl::read_Triangle_node(filename, mesh_nodes_, gctl::Packed, &node_markers_);
    gctl::read_Triangle_element(filename, mesh_elems_, mesh_nodes_, gctl::Packed, &tmp_markers);

    elem_num_ = mesh_elems_.size();
    node_num_ = mesh_nodes_.size();

    elem_markers_.resize(elem_num_);
    for (int i = 0; i < elem_num_; ++i)
    {
        elem_markers_[i] = (int) tmp_markers[i][0];
    }
    tmp_markers.clear();

    std::clog << "done\n";
    std::clog << "triangle number = " << elem_num_ << ", node number = " << node_num_ << std::endl;
    return;
}

// Output a Gmsh file
void mt2d::write_model(std::string filename, bool debug)
{
    std::clog << "Writing Gmsh file ... ";
    meshio_.init_file(filename, gctl::Output);
    meshio_.set_packed(gctl::NotPacked, gctl::Output);
    meshio_.save_mesh(mesh_elems_, mesh_nodes_);
    meshio_.save_data("vertex marker", node_markers_, gctl::NodeData);
    meshio_.save_data("element marker", elem_markers_, gctl::ElemData);

    if (debug)
    {
        gctl::array<double> d1(node_num_, 0.0), d2(node_num_, 0.0);
        for (int i = 0; i < node_num_; ++i)
        {
            d1[i] = x[i].real();
            d2[i] = x[i].imag();
        }

        if (work_mode_ == TE)
        {
            meshio_.save_data("Ex_real", d1, gctl::NodeData);
            meshio_.save_data("Ex_imag", d2, gctl::NodeData);

            gctl::array<double> ta(node_num_, 0.0);
            gctl::array<double> d3(node_num_, 0.0), d4(node_num_, 0.0);

            for (int i = 0; i < elem_num_; ++i)
            {
                for (int j = 0; j < 3; ++j)
                {
                    ta[mesh_elems_[i].vert[j]->id] += mesh_elems_[i].area();
                }
            }

            complex_d Hy;
            for (int i = 0; i < elem_num_; ++i)
            {
                for (int j = 0; j < 3; ++j)
                {
                    Hy = _zero;
                    for (size_t k = 0; k < 3; k++)
                    {
                        Hy += x[mesh_elems_[i].vert[k]->id]*lsf_.triangle(mesh_elems_[i].vert[j]->x, mesh_elems_[i].vert[j]->y, 
                            mesh_elems_[i].vert[0]->x, mesh_elems_[i].vert[1]->x, mesh_elems_[i].vert[2]->x, 
                            mesh_elems_[i].vert[0]->y, mesh_elems_[i].vert[1]->y, mesh_elems_[i].vert[2]->y, k, nfem_sf::Dy);
                    }
                    Hy /= (_1i*omega_*_mu0);

                    d3[mesh_elems_[i].vert[j]->id] += Hy.real()*mesh_elems_[i].area()/ta[mesh_elems_[i].vert[j]->id];
                    d4[mesh_elems_[i].vert[j]->id] += Hy.imag()*mesh_elems_[i].area()/ta[mesh_elems_[i].vert[j]->id];
                }
            }

            meshio_.save_data("Hy_real", d3, gctl::NodeData);
            meshio_.save_data("Hy_imag", d4, gctl::NodeData);
        }
        else // TM
        {
            meshio_.save_data("Hx_real", d1, gctl::NodeData);
            meshio_.save_data("Hx_imag", d2, gctl::NodeData);

            gctl::array<double> ta(node_num_, 0.0);
            gctl::array<double> d3(node_num_, 0.0), d4(node_num_, 0.0);

            for (int i = 0; i < elem_num_; ++i)
            {
                for (int j = 0; j < 3; ++j)
                {
                    ta[mesh_elems_[i].vert[j]->id] += mesh_elems_[i].area();
                }
            }

            complex_d Ey;
            double elem_cndt;
            for (int i = 0; i < elem_num_; ++i)
            {
                get_element_conductivity(elem_markers_[i], elem_cndt);

                for (int j = 0; j < 3; ++j)
                {
                    Ey = _zero;
                    for (size_t k = 0; k < 3; k++)
                    {
                        Ey += -1.0*(1.0/elem_cndt)*x[mesh_elems_[i].vert[k]->id]*lsf_.triangle(mesh_elems_[i].vert[j]->x, mesh_elems_[i].vert[j]->y, 
                            mesh_elems_[i].vert[0]->x, mesh_elems_[i].vert[1]->x, mesh_elems_[i].vert[2]->x, 
                            mesh_elems_[i].vert[0]->y, mesh_elems_[i].vert[1]->y, mesh_elems_[i].vert[2]->y, k, nfem_sf::Dy);
                    }

                    d3[mesh_elems_[i].vert[j]->id] += Ey.real()*mesh_elems_[i].area()/ta[mesh_elems_[i].vert[j]->id];
                    d4[mesh_elems_[i].vert[j]->id] += Ey.imag()*mesh_elems_[i].area()/ta[mesh_elems_[i].vert[j]->id];
                }
            }

            meshio_.save_data("Ey_real", d3, gctl::NodeData);
            meshio_.save_data("Ey_imag", d4, gctl::NodeData);
        }
    }
    
    std::clog << "done\n";
    return;
}

void mt2d::write_observation(std::string filename)
{
    std::clog << "Writing Observation file ... ";
    std::ofstream outfile;
    gctl::open_outfile(outfile, filename, ".obs.txt");
    for (size_t i = 0; i < obs_num_; i++)
    {
        outfile << obs_site_[i].vert->x << " " << obs_site_[i].vert->y << " " 
            << obs_site_[i].app_resist << " " << obs_site_[i].phase << std::endl;
    }
    outfile.close();
    std::clog << "done\n";
    return;
}

void mt2d::get_element_conductivity(int tag_num, double &conduct, bool ref_cndt)
{
    for (int i = 0; i < secondary_model_.size(); ++i) // 数量有限直接搜索即可
    {
        if (secondary_model_[i].tag == tag_num)
        {
            if (ref_cndt) conduct = secondary_model_[i].ref_cndt;
            else conduct = secondary_model_[i].cndt;
            return;
        }
    }

    throw gctl::runtime_error("No secondary model found by the tag number: "+std::to_string(tag_num));
    return;
}

void mt2d::assemble_linear_system()
{
    std::clog << "Assembling the linear system ... ";

    std::vector<triplt_cd> kernel_entries;
    kernel_entries.reserve(node_num_*3);

    g2d_tri_.initiate(2); // 初始化高斯积分阶数

    // 计算线性方程组的左端项
    A.resize(node_num_, node_num_);
    A.setZero();

    unsigned int uv_id[2];
    double elem_cndt, aij_p1, aij_p2;
    complex_d tmp_entry;
    gctl::vertex2dc* tmp_vert[3];
    // Loop all elements
    for (int e = 0; e < elem_num_; ++e)
    {
        tmp_vert[0] = mesh_elems_[e].vert[0];
        tmp_vert[1] = mesh_elems_[e].vert[1];
        tmp_vert[2] = mesh_elems_[e].vert[2];

        get_element_conductivity(elem_markers_[e], elem_cndt);

        for (int i = 0; i < 3; i++)
        {
            for (int j = 0; j < 3; j++)
            {
                uv_id[0] = i; uv_id[1] = j;

                aij_p1 = g2d_tri_.integral(tmp_vert[0]->x, tmp_vert[1]->x, tmp_vert[2]->x, 
                    tmp_vert[0]->y, tmp_vert[1]->y, tmp_vert[2]->y, _Aij_part1, &uv_id[0], this);

                aij_p2 = g2d_tri_.integral(tmp_vert[0]->x, tmp_vert[1]->x, tmp_vert[2]->x, 
                    tmp_vert[0]->y, tmp_vert[1]->y, tmp_vert[2]->y, _Aij_part2, &uv_id[0], this);

                if (work_mode_ == TE)
                {
                    tmp_entry.real(aij_p1);
                    tmp_entry.imag(-1.0*omega_*_mu0*elem_cndt*aij_p2);
                }
                else // TM
                {
                    tmp_entry.real(aij_p1/elem_cndt);
                    tmp_entry.imag(-1.0*omega_*_mu0*aij_p2);
                }

                kernel_entries.push_back(triplt_cd(tmp_vert[i]->id, tmp_vert[j]->id, tmp_entry));
            }
        }
    }

    A.setFromTriplets(kernel_entries.begin(), kernel_entries.end());
    kernel_entries.clear();

    x.resize(node_num_);
    x.setZero();

    // 计算线性方程组的右端项 对于一次场方法右端项为0
    B.resize(node_num_);
    B.setZero();

    // 计算线性方程组的右端项 仅在使用二次场方法时计算
    unsigned int u_id;
    double elem_cndt_ref;
    complex_d c_factor;
    for (int e = 0; e < elem_num_; ++e)
    {
        tmp_vert[0] = mesh_elems_[e].vert[0];
        tmp_vert[1] = mesh_elems_[e].vert[1];
        tmp_vert[2] = mesh_elems_[e].vert[2];

        get_element_conductivity(elem_markers_[e], elem_cndt);
        get_element_conductivity(elem_markers_[e], elem_cndt_ref, true);

        if (elem_cndt != elem_cndt_ref)
        {
            if (work_mode_ == TE) c_factor = _1i*omega_*_mu0*(elem_cndt - elem_cndt_ref);
            else c_factor = (1.0 - elem_cndt_ref/elem_cndt);

            for (int i = 0; i < 3; i++)
            {
                u_id = i;
                B[mesh_elems_[e].vert[i]->id] += c_factor*g2d_tri_.integral(tmp_vert[0]->x, tmp_vert[1]->x, tmp_vert[2]->x, 
                    tmp_vert[0]->y, tmp_vert[1]->y, tmp_vert[2]->y, _Bi, &u_id, this);
            }
        }
    }

    std::clog << "done\n";
    return;
}

void mt2d::apply_boundary_condition()
{
    std::clog << "Applying boundary condition ... ";

    for (int i = 0; i < node_num_; ++i)
    {
        if (node_markers_[i] == bound_marker_) // 边界上的顶点
        {
            B.coeffRef(i) = _zero; // 二次场方法计算中假设边界上的二次场已经衰减为0
        }
    }

    // 采用划0置1的方法处理第一类边界条件 （稀疏矩阵处理起来很快）

    for (int i = 0; i < node_num_; ++i) // 外层是列循环
    {
        if (node_markers_[i] == bound_marker_)
        {
            for (Eigen::SparseMatrix<complex_d>::InnerIterator it(A, i); it; ++it) // 只循环非零值 行循环
            {
                if (node_markers_[it.row()] != bound_marker_)  // 非边界上的顶点
                {
                    B[it.row()] -= it.value() * B[i];
                }
            }
        }
    }

    for (int i = 0; i < node_num_; i++) // 处理行
    {
        for (Eigen::SparseMatrix<complex_d>::InnerIterator it(A, i); it; ++it)
        {
            if (node_markers_[it.row()] == bound_marker_)
            {
                it.valueRef() = _zero;
            }
        }
    }

    for (int i = 0; i < node_num_; ++i) // 处理列
    {
        if (node_markers_[i] == bound_marker_)
        {
            for (Eigen::SparseMatrix<complex_d>::InnerIterator it(A, i); it; ++it)
            {
                it.valueRef() = _zero;
                if (it.row() == it.col()) it.valueRef() = _one; // 置1
            }
        }
    }

    A.makeCompressed();
    // 边界条件处理完成

    std::clog << "done\n";
    return;
}

void mt2d::solve()
{
    std::clog << "Solving the PDE ... ";
    // Call SparseLU to solve the linear system
    Eigen::SparseLU<spmat_cd> solver;
    solver.analyzePattern(A);
    solver.factorize(A);
    x = solver.solve(B);
    std::clog << "done\n";
    return;
}

void mt2d::post_process()
{
    // add primary field
    if (work_mode_ == TE)
    {
        for (int i = 0; i < node_num_; ++i)
        {
            le_.calculate_EH(-1.0*mesh_nodes_[i].y, mesh_nodes_[i].x);
            x[i] += le_.get_Ex();
        }
    }
    else // TM
    {
        for (int i = 0; i < node_num_; ++i)
        {
            le_.calculate_EH(-1.0*mesh_nodes_[i].y, mesh_nodes_[i].x);
            x[i] += le_.get_Hx();
        }
    }

    // Sort out MT stations
    mt_station tmp_station;
    for (size_t i = 0; i < node_num_; i++)
    {
        if (node_markers_[i] == station_marker_)
        {
            tmp_station.vert = mesh_nodes_.get(i);
            tmp_station.app_resist = 0.0;
            tmp_station.phase = 0.0;

            obs_site_.push_back(tmp_station);
        }
    }

    obs_num_ = obs_site_.size();

    // Find triangles that include mt stations
    std::vector<mt_station>::iterator stat_iter;
    for (size_t i = 0; i < elem_num_; i++)
    {
        if (elem_markers_[i] != air_marker_) // exclude the air layer to obtain smoother results
        {
            for (size_t j = 0; j < 3; j++)
            {
                if (node_markers_[mesh_elems_[i].vert[j]->id] == station_marker_)
                {
                    for (size_t k = 0; k < obs_num_; k++)
                    {
                        if (obs_site_[k].vert == mesh_elems_[i].vert[j])
                        {
                            obs_site_[k].host_tris.push_back(mesh_elems_.get(i));
                            break;
                        }
                    }
                }
            }
        }
    }

    // Calculate observations
    double total_area, elem_cndt;
    complex_d Ey, Hy;

    if (work_mode_ == TE)
    {
        for (size_t i = 0; i < obs_num_; i++)
        {
            total_area = 0.0;
            for (size_t j = 0; j < obs_site_[i].host_tris.size(); j++)
            {
                total_area += obs_site_[i].host_tris[j]->area();
            }

            Hy = _zero;
            for (size_t j = 0; j < obs_site_[i].host_tris.size(); j++)
            {
                // 下面其实就是一个平面函数插值过程
                for (size_t k = 0; k < 3; k++)
                {
                    Hy += (obs_site_[i].host_tris[j]->area()/total_area)*x[obs_site_[i].host_tris[j]->vert[k]->id]
                        *lsf_.triangle(obs_site_[i].vert->x, obs_site_[i].vert->y, 
                        obs_site_[i].host_tris[j]->vert[0]->x, obs_site_[i].host_tris[j]->vert[1]->x, obs_site_[i].host_tris[j]->vert[2]->x, 
                        obs_site_[i].host_tris[j]->vert[0]->y, obs_site_[i].host_tris[j]->vert[1]->y, obs_site_[i].host_tris[j]->vert[2]->y, 
                        k, nfem_sf::Dy);
                }
            }
            Hy /= -1.0*_1i*omega_*_mu0; // 注意这里纵轴向上与层状模型坐标轴方向相反 Dy方向导数应该乘-1

            obs_site_[i].app_resist = std::norm(x[obs_site_[i].vert->id]/Hy)/(omega_*_mu0);
            obs_site_[i].phase = std::arg(x[obs_site_[i].vert->id]/Hy)*180.0/pi;
        }
    }
    else // TM
    {
        for (size_t i = 0; i < obs_num_; i++)
        {
            total_area = 0.0;
            for (size_t j = 0; j < obs_site_[i].host_tris.size(); j++)
            {
                total_area += obs_site_[i].host_tris[j]->area();
            }

            Ey = _zero;
            for (size_t j = 0; j < obs_site_[i].host_tris.size(); j++)
            {
                get_element_conductivity(elem_markers_[obs_site_[i].host_tris[j]->id], elem_cndt);
                // 下面其实就是一个平面函数插值过程
                for (size_t k = 0; k < 3; k++)
                {
                    Ey += -1.0*(1.0/elem_cndt)*(obs_site_[i].host_tris[j]->area()/total_area)*x[obs_site_[i].host_tris[j]->vert[k]->id]
                        *lsf_.triangle(obs_site_[i].vert->x, obs_site_[i].vert->y, 
                        obs_site_[i].host_tris[j]->vert[0]->x, obs_site_[i].host_tris[j]->vert[1]->x, obs_site_[i].host_tris[j]->vert[2]->x, 
                        obs_site_[i].host_tris[j]->vert[0]->y, obs_site_[i].host_tris[j]->vert[1]->y, obs_site_[i].host_tris[j]->vert[2]->y, 
                        k, nfem_sf::Dy); // 注意这里纵轴向上与层状模型坐标轴方向相反 Dy方向导数应该乘-1
                }
            }

            obs_site_[i].app_resist = std::norm(Ey/x[obs_site_[i].vert->id])/(omega_*_mu0);
            obs_site_[i].phase = std::arg(Ey/x[obs_site_[i].vert->id])*180.0/pi;
        }
    }
    return;
}

void mt2d::run(int argc, char *argv[])
{
    // Read configuration file
    read_config(argv[1]);
    // Read model file
    read_model(tri_file_);

    if (work_mode_ != Modeling)
    {
        // Prepare the linear system
        assemble_linear_system();
        // Apply boundary condition
        apply_boundary_condition();
        // Solve the linear system
        solve();
        // Calculate observations
        post_process();
        // Save Gmsh model file
        write_model(tri_file_);
        // Write observations
        write_observation(tri_file_);
    }
    else
    {
        // Save Gmsh model file
        write_model(tri_file_, false);
    }
    return;
}

int main(int argc, char *argv[]) try
{
    if (argc < 2)
    {
        std::clog << "Usage: ./mt2d <config_file>\n";
        return -1;
    }

    mt2d instance;
    instance.run(argc, argv);
    return 0;
}
catch (std::exception &e)
{
    GCTL_ShowWhatError(e.what(), GCTL_ERROR_ERROR, 0, 0, 0);
}
